/*******************************************************************************
Paper title: Does early nutrition predict cognitive skills during later childhood?
Evidence from two developing countries
Alan Sanchez, Marta Favara, Margaret Sheridan, Jere Behrman.
Created: 8 OCt 2020
This version: 4 Dec 2023 
Content: Figure 1, Figure B1, Figure B2, Figure B4
*******************************************************************************/
			
clear all
local user Alan

if "`user'"=="Alan" {
global outdata	C:\Users\alans\Dropbox\_NIHR21_Proposal\NIH Papers\Paper 1_nutrition&EF\paper\WD\Replication files\Data
global output   C:\Users\alans\Dropbox\_NIHR21_Proposal\NIH Papers\Paper 1_nutrition&EF\paper\WD\Replication files\Output
}

global control0   chage_r4 chage_r4_2 female dm_educ2 dm_educ3 dm_educ4 urban_c hhsize_c
global country1   mom_spanish
global country2   mom_oromifa mom_tigrina mom_other
global control2q  wi_cq2 wi_cq3 wi_cq4 wi_cq5
global alltask 	  hr2 hr3 wkend /*practices*/
global cluster_r1 sib_y sib_o /*dclustid_r11- dclustid_r120*/

global control0xf  i.female*chage_r4 i.female*chage_r4_2 i.female*dm_educ2 ///
                   i.female*dm_educ3 i.female*dm_educ4 i.female*urban_c i.female*hhsize_c
global control2qf  i.female*wi_cq2 i.female*wi_cq3 i.female*wi_cq4 i.female*wi_cq5  
global country1f   i.female*mom_spanish 
global country2f   i.female*mom_oromifa i.female*mom_tigrina i.female*mom_other
global alltaskf    i.female*hr2 i.female*hr3 i.female*wkend /*i.female*practices*/
global cluster_r1f i.female*sib_y i.female*sib_o

global control0xa  chage_r4xdm_educ2 chage_r4xdm_educ3 chage_r4xdm_educ4 ///
                   chage_r4xurban_c chage_r4xhhsize_c
global control2qa  chage_r4xwi_cq2 chage_r4xwi_cq3 chage_r4xwi_cq4 chage_r4xwi_cq5  
global country1a   chage_r4xmom_spanish 
global country2a   chage_r4xmom_oromifa chage_r4xmom_tigrina chage_r4xmom_other
global alltaska    chage_r4xhr2 chage_r4xhr3 chage_r4xwkend /*chage_r4xpractices*/
global cluster_r1a chage_r4xsib_y chage_r4xsib_o

global task1 	task1
global task2 	itask2_a
global task3 	itask3
global task4 	itask4	

global btask1 	btask1
global btask2 	bitask2_a
global btask3 	bitask3
global btask4 	bitask4	

global sample1  zhfa_c<=6 & zhfa_c>=-6
global sample2  zhfa_c<=6 & zhfa_c>=-6

global absorb   clustid_r1
global indlevel pride_index_r4 agency_index_r4 sesteem_index_r4 ppvtz_r4 missing_ppvt

*---------------------------------------------------------------------------------------------------------------------*---------------------------------------------------------------------------------------------------------------------*------------------------------------------------------------------------------------------------------------------*/
use "$outdata\pe_et_childlevel_230821", clear

bysort childid sibling: egen n=max(task) /* keep only paired-siblings who answer all tasks*/
drop if n!=4 /*22 obs*/
drop n
est clear

gen sib_y=0
gen sib_o=0
replace sib_y=1 if sibling==1 & sibyounger==1
replace sib_o=1 if sibling==1 & sibyounger==0 

gen zhfa_c=.
replace zhfa_c=zhfa_r3 if sibling==1
replace zhfa_c=zhfa_r2 if sibling==0

gen stunting_c=.
replace stunting_c=stunting_r3 if sibling==1
replace stunting_c=stunting_r2 if sibling==0

gen zhfaxsib_y=zhfa_c*sib_y
gen zhfaxsib_o=zhfa_c*sib_o
gen zhfaxage=zhfa_c*chage_r4

gen stuntingxsib_y=stunting_c*sib_y
gen stuntingxsib_o=stunting_c*sib_o
gen stuntingxage  =stunting_c*chage_r4

gen zhfaxfem=zhfa_c*female

gen stuntingxfem=stunting_c*female

gen chage_r4_2=chage_r4*chage_r4

label var stunting_c     "Stunted"
label var stuntingxsib_y "Stunted x younger sib"
label var stuntingxsib_o "Stunted x older sib"
label var stuntingxfem   "Stunted x female"

label var zhfa_c         "Height for age"
label var zhfaxsib_y     "Height for age x younger sib"
label var zhfaxsib_o     "Height for age x older sib"
label var zhfaxfem       "Height for age x female"

label var chage_r4     "Age in months, r4"
label var chage_r4_2   "Age in months squared, r4"
label var female       "Child is female"
label var dm_educ2     "Maternal edu: complete primary"
label var dm_educ3     "Maternal edu: complete secondary"
label var dm_educ4     "Maternal edu: complete tertiary"
label var urban_r1     "Urban area, r1"
label var hhsize_r2    "Household size, r2" 
label var mom_spanish  "Maternal native tongue: spanish"
label var mom_oromifa  "Maternal native tongue: oromifah" 
label var mom_tigrina  "Maternal native tongue: tigrina"
label var mom_other    "Maternal native tongue: other"
label var wi_r1q2      "Wealth index quintile 2"
label var wi_r1q3      "Wealth index quintile 3"
label var wi_r1q4      "Wealth index quintile 4"
label var wi_r1q5      "Wealth index quintile 5"

*****************************************************
*FIGURE B1 (Appendix)
*****************************************************
twoway (hist chage_r4 if sibling==0 & country==1, percent bfcolor(none) blcolor(orange) width(2)) ///
(hist chage_r4 if sibyounger==1 & country==1, percent bfcolor(none) blcolor(navy) width(2)) ///
	   , xline(10) title(Young Lives children in Peru) xtitle(Age in months)  ///
legend(label(1 "Index children") label(2 "Younger siblings") style(none)) 
graph save peru_haz.gph, replace

twoway (hist chage_r4 if sibling==0 & country==2, percent bfcolor(none) blcolor(orange) width(2)) ///
(hist chage_r4 if sibyounger==1 & country==2, percent bfcolor(none) blcolor(navy) width(2)) ///
(hist chage_r4 if sibyounger==0 & country==2, percent bfcolor(none) blcolor(pink) width(2)) ///
	   , xline(10) title(Young Lives children in Ethiopia) xtitle(Age in months)  ///
legend(label(1 "Index children") label(2 "Younger siblings") label(3 "Older siblings") style(none)) 
graph save ethiopia_haz.gph, replace
graph combine peru_haz.gph ethiopia_haz.gph, col(1) xcommon ycommon 
graph save "$output\figB1.gph", replace
*END OF FIGURE B1

*****************************************************
*FIGURE B2 (Appendix) 
*****************************************************
*Average HAZ by age in semesters 

preserve
duplicates drop childid sibling, force
keep country childid sibling sibyounger chage_r1 chage_r2 chage_r3 chage_r4 chage_r5 zhfa_r1 zhfa_r2 zhfa_r3 zhfa_r4 zhfa_r5 
reshape long chage_r zhfa_r chsex, i(childid sibling) j(round)
drop if zhfa_r<=-6
drop if zhfa_r>=6
gen sem=.
replace sem=1 if chage_r>=6   & chage_r<=12
replace sem=2 if chage_r>=13  & chage_r<=18
replace sem=3 if chage_r>=19  & chage_r<=24
replace sem=4 if chage_r>=25  & chage_r<=30

replace sem=5 if chage_r>=31  & chage_r<=36
replace sem=6 if chage_r>=37  & chage_r<=42
replace sem=7 if chage_r>=43  & chage_r<=48
replace sem=8 if chage_r>=49  & chage_r<=54
replace sem=9 if chage_r>=55  & chage_r<=60

replace sem=10 if chage_r>=61  & chage_r<=66
replace sem=11 if chage_r>=67  & chage_r<=72
replace sem=12 if chage_r>=73  & chage_r<=78
replace sem=13 if chage_r>=79  & chage_r<=84
replace sem=14 if chage_r>=85  & chage_r<=90

replace sem=15 if chage_r>=91  & chage_r<=96
replace sem=16 if chage_r>=97  & chage_r<=102
replace sem=17 if chage_r>=103 & chage_r<=108
replace sem=18 if chage_r>=109 & chage_r<=114
replace sem=19 if chage_r>=115 & chage_r<=120

replace sem=20 if chage_r>=121 & chage_r<=126
replace sem=21 if chage_r>=127 & chage_r<=132
replace sem=22 if chage_r>=133 & chage_r<=138
replace sem=23 if chage_r>=139 & chage_r<=144
replace sem=24 if chage_r>=145 & chage_r<=150

replace sem=25 if chage_r>=151 & chage_r<=156
replace sem=26 if chage_r>=157 & chage_r<=162
replace sem=27 if chage_r>=163 & chage_r<=168
replace sem=28 if chage_r>=169 & chage_r<=174
replace sem=29 if chage_r>=175 & chage_r<=180

replace sem=30 if chage_r>=181 & chage_r<=186
replace sem=31 if chage_r>=187 & chage_r<=192
replace sem=32 if chage_r>=193 & chage_r<=198
replace sem=33 if chage_r>=199 & chage_r<=204
replace sem=34 if chage_r>=205 & chage_r<=210

replace sem=35 if chage_r>=211 & chage_r<=216
replace sem=36 if chage_r>=217 & chage_r<=222
replace sem=37 if chage_r>=223 & chage_r<=228
replace sem=38 if chage_r>=229 & chage_r<=234
replace sem=39 if chage_r>=235 & chage_r<=240

drop if round==1

collapse (mean) meanhaz=zhfa_r (count) counthaz=zhfa_r, by(country sem sibling)
drop if counthaz<10

twoway (scatter meanhaz sem if sibling==0 & country==1, mcolor(none) mlcolor(green) mfcolor(green)) ///
       (scatter meanhaz sem if sibling==1 & country==1, mcolor(none) mlcolor(black) mfcolor(white)) ///
, xline(10) title(Young Lives children in Peru) xtitle(Age in semesters) ///
ytitle(Average HAZ) ///
legend(label(1 "Index children") label(2 "Siblings")) 
graph save peru_haz.gph, replace

twoway (scatter meanhaz sem if sibling==0 & country==2, mcolor(none) mlcolor(green) mfcolor(green)) ///
       (scatter meanhaz sem if sibling==1 & country==2, mcolor(none) mlcolor(black) mfcolor(white)) ///
, xline(10) title(Young Lives children in Ethiopia) xtitle(Age in semesters) ///
ytitle(Average HAZ) ///
legend(label(1 "Index children") label(2 "Siblings")) 
graph save ethiopia_haz.gph, replace

graph combine peru_haz.gph ethiopia_haz.gph, col(1) xcommon ycommon
graph save "$output\figB2.gph", replace

restore
*END OF FIGURE B2

***Index children
gen     zhfa_5 =zhfa_r2     if (zhfa_r2>=-6 & zhfa_r2<=6) & sibling==0                 & country==2 
gen     zhfa_12=zhfa_r4     if (zhfa_r4>=-6 & zhfa_r4<=6) & sibling==0                 & country==2 

***Older siblings
replace zhfa_12=zhfa_r3 if (zhfa_r3>=-6 | zhfa_r3<=6) & sibling==1 & sibyounger==0 & country==2 

***PPVT
gen missing_ppvt=.
replace missing_ppvt=0 if country==2 & ppvtz_r4!=.
replace missing_ppvt=1 if country==2 & ppvtz_r4==. 
label var missing_ppvt "Missing PPVT"
tab missing_ppvt
replace ppvtz_r4=10000 if ppvtz_r4==.

***Socio emotional outcomes
gen missing_pride_r4=.
replace missing_pride=0 if country==2 & pride_index_r4!=.
replace missing_pride=1 if country==2 & pride_index_r4==. 
label var missing_pride "Missing pride index"
tab missing_pride
replace pride_index_r4=10000 if pride_index_r4==.

gen missing_agency_r4=.
replace missing_agency=0 if country==2 & agency_index_r4!=.
replace missing_agency=1 if country==2 & agency_index_r4==. 
label var missing_agency "Missing agency index"
tab missing_agency
replace agency_index_r4=10000 if agency_index_r4==.

gen missing_sesteem_r4=.
replace missing_sesteem=0 if country==2 & sesteem_index_r4!=.
replace missing_sesteem=1 if country==2 & sesteem_index_r4==. 
label var missing_sesteem "Missing self-esteem index"
tab missing_sesteem
replace sesteem_index_r4=10000 if sesteem_index_r4==.

foreach x of varlist zhfa_5 zhfa_12 pride_index_r4 agency_index_r4 sesteem_index_r4 ppvtz_r4 missing_ppvt ///
hsleep_r4 hcare_r4 hchore_r4 htask_r4 hwork_r4 hschool_r4 hstudy_r4 hplay_r4 female ///
dm_educ2 dm_educ3 dm_educ4 urban_c hhsize_c mom_oromifa mom_tigrina mom_other wi_cq2 wi_cq3 wi_cq4 wi_cq5 {
	gen `x'xage  =. 
	gen `x'xage2 =.
	replace `x'xage  =`x'*chage_r2           if sibling==0 
	replace `x'xage  =`x'*chage_r3           if sibling==1
	replace `x'xage2 =`x'*chage_r2*chage_r2  if sibling==0
	replace `x'xage2 =`x'*chage_r3*chage_r3  if sibling==1
}

forvalues x= 1/20 {
	gen clust`x'_=0
	replace clust`x'_=1 if clustid_r1==`x'
}

********************************************************************************
*** Adjusted HAZ for older siblings in Ethiopia
********************************************************************************
gen zhfa_or=zhfa_c
gen stunting_o=stunting_c
sort  country childid sibling

merge country childid sibling using "$outdata\lasso_prediction_older.dta"
drop _merge

replace zhfa_c      =zhfa_5_p               if sibling==1 & sibyounger==0 & country==2
replace zhfaxsib_y  =zhfa_c*sib_y           if sibling==1 & sibyounger==0 & country==2
replace zhfaxsib_o  =zhfa_c*sib_o           if sibling==1 & sibyounger==0 & country==2
replace zhfaxage    =zhfa_c*chage_r4        if sibling==1 & sibyounger==0 & country==2

tab     stunting_c                 		    if sibling==1 & sibyounger==0 & country==2 & task1!=.
replace stunting_c    =1          		    if zhfa_c<-2  & zhfa_c!=. & sibling==1 & sibyounger==0 & country==2
replace stunting_c    =0          		    if zhfa_c>=-2 & zhfa_c!=. & sibling==1 & sibyounger==0 & country==2
replace stuntingxsib_y=stunting_c*sib_y  	if sibling==1 & sibyounger==0 & country==2
replace stuntingxsib_o=stunting_c*sib_o   	if sibling==1 & sibyounger==0 & country==2
replace stuntingxage  =stunting_c*chage_r4 	if sibling==1 & sibyounger==0 & country==2
tab     stunting_c                   		if sibling==1 & sibyounger==0 & country==2 & task1!=.

********************************************************************************
*** Adjusted for index children and younger siblings in Peru & Ethiopia
********************************************************************************
sort  country childid sibling
merge country childid sibling using "$outdata\prediction_index_younger.dta"
drop _merge

replace zhfa_c      =zhfa_5_p1              if sibyounger!=0
replace zhfaxsib_y  =zhfa_c*sib_y           if sibyounger!=0
replace zhfaxsib_o  =zhfa_c*sib_o           if sibyounger!=0
replace zhfaxage    =zhfa_c*chage_r4        if sibyounger!=0

tab     stunting_c country         		    if sibyounger!=0 & task1!=.
replace stunting_c    =1          		    if zhfa_c<-2  & zhfa_c!=. & sibyounger!=0
replace stunting_c    =0          		    if zhfa_c>=-2 & zhfa_c!=. & sibyounger!=0
replace stuntingxsib_y=stunting_c*sib_y  	if sibyounger!=0
replace stuntingxsib_o=stunting_c*sib_o   	if sibyounger!=0
replace stuntingxage  =stunting_c*chage_r4 	if sibyounger!=0
tab     stunting_c                   		if sibyounger!=0 & task1!=.

********************************************************************
***FIGURE 1
********************************************************************
twoway (hist zhfa_c if country==1 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(orange)) ///
(hist zhfa_o if country==1 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(navy) title("Young Lives pooled sample in Peru") ///
xline(-2) legend(lab(1 "Adjusted HAZ") lab(2 "Observed HAZ")))
graph save peru_haz.gph, replace

twoway (hist zhfa_c if country==2 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(orange)) ///
(hist zhfa_o if country==2 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(navy) title("Young Lives pooled sample in Ethiopia") ///
xline(-2) legend(lab(1 "Adjusted HAZ") lab(2 "Observed HAZ")))
graph save ethiopia_haz.gph, replace

graph combine peru_haz.gph ethiopia_haz.gph, col(1) xcommon ycommon
graph save "$output\fig1.gph", replace
*END OF FIGURE 1

*****************************************************
*FIGURE B4, part I
*****************************************************
twoway (hist zhfa_c if country==1 & sibling==0 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(orange)) ///
(hist zhfa_o if country==1 & sibling==0 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(navy) title("Predicted versus observed: Peru (index children)") ///
xline(-2) legend(lab(1 "Adjusted HAZ") lab(2 "Observed HAZ")))
graph save peru_haz.gph, replace

twoway (hist zhfa_c if country==2 & sibling==0 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(orange)) ///
(hist zhfa_o if country==2 & sibling==0 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(navy) title("Predicted versus observed: Ethiopia (index children)") ///
xline(-2) legend(lab(1 "Adjusted HAZ") lab(2 "Observed HAZ")))
graph save ethiopia_haz.gph, replace

graph combine peru_haz.gph ethiopia_haz.gph, col(1) xcommon ycommon
graph save "$output\figB4I.gph", replace

*****************************************************
*FIGURE B4, part II
*****************************************************
twoway (hist zhfa_c if country==1 & sibyounger==1 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(orange)) ///
(hist zhfa_o if country==1 & sibyounger==1 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(navy) title("Predicted versus observed: Peru (younger siblings)") ///
xline(-2) legend(lab(1 "Adjusted HAZ") lab(2 "Observed HAZ")))
graph save peru_haz.gph, replace

twoway (hist zhfa_c if country==2 & sibyounger==1 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(orange)) ///
(hist zhfa_o if country==2 & sibyounger==1 & $sample1 & task1!=., ///
percent width(0.04) bfcolor(none) blcolor(navy) title("Predicted versus observed: Ethiopia (younger siblings)") ///
xline(-2) legend(lab(1 "Adjusted HAZ") lab(2 "Observed HAZ")))
graph save ethiopia_haz.gph, replace

graph combine peru_haz.gph ethiopia_haz.gph, col(1) xcommon ycommon
graph save "$output\figB4II.gph", replace
*END OF FIGURE B4